Heterogeneous animal group models and their group-level alignment dynamics; 

an equation- free approach 
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We study coarse-grained (group- level) alignment dynamics of individual-based animal group mod- 
els for heterogeneous populations consisting of informed (on preferred directions) and uninformed 
individuals. The orientation of each individual is characterized by an angle, whose dynamics are 
nonlinearly coupled with those of all the other individuals, with an explicit dependence on the 
difference between the individual's orientation and the instantaneous average direction. Choosing 
convenient coarse-grained variables (suggested by uncertainty quantification methods) that account 
for rapidly developing correlations during initial transients, we perform efficient computations of 
coarse-grained steady states and their bifurcation analysis. We circumvent the derivation of coarse- 
grained governing equations, following an equation-free computational approach. 
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I. INTRODUCTION 

Coordinated motions and pattern formation have been 
studied for a wide range of biological organisms, from 
bacteria and amoebae to fish, from birds and wildebeest 
to humans (Ben- Jacob et al., 2000; Deneubourg et al., 
2001; Parrishet al., 2002; Partridge, 1982; Wilson, 1975). 
Animal groups often behave as if they have a single mind, 
displaying remarkable self-organized behavior. At one 
extreme, the individuals seem to need little information 
transfer (e.g. fish schools), while at the other end the 
information exchange occurs in highly integrated ways 
through long-term associations among the individuals 
(e.g. honeybee hives and human communities). Con- 
trolling such an organized behavior in groups of artifi- 
cial objects, including autonomous underwater vehicles 
(Leonard et al., 2006) and groups of autonomous agents 
(Jadbabaie et al., 2003), has received extensive attention 
in contemporary control theory. Challenges in both natu- 
ral and engineering settings involve understanding which 
patterns emerge from the interaction among individual 
agents. 

Select laboratory experiments have shed some light on 
the schooling mechanism (Hunter, 1966; van 01st and 
Hunter, 1970; Partridge and Pitcher, 1979, 1980; Pitcher 
et al., 1976). It still remains unclear, however, how 
the individual- level behavior and group- level ("macro- 
scopic", or coarse-grained) patterns are related. More 
precise experiments using three-dimensional tracking of 
every individual in a population should lead to better 
understanding of this linkage. An ultimate experimental 
study with precise control of every relevant detail may 
not be possible, yet appropriate mathematical models 
would provide a venue to establish behavioral cause, as 
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one can consider different hypothetical individual-level 
interaction rules selectively (see e.g., Flierl et al., 1999). 

Several different individual-based models have been 
proposed, which reproduce certain types of collective be- 
havior in animal groups (e.g., see Aoki, 1982; Reynolds, 
1987; Deneubourg and Goss, 1989). Self-organization 
emerges also in a wide spectrum of physical and chemi- 
cal systems, some of which (e.g., crystals and ferromag- 
netic materials) exhibit apparent similarities with emer- 
gent patterns observed in animal groups. Vicsek et al. 
(1995) have introduced a discrete-time model of self- 
driven particles, or self-propelled particles (SPP), based 
on near-neighbor rules that are similar with those in 
the ferromagnetic XY model (Kosterlitz and Thouless, 
1973). The authors analyzed statistical properties of the 
model, including phase transition and scaling (Vicsek et 
al., 1995). A long-range interaction has been incorpo- 
rated into the SPP model (Mikhailov and Zanette, 1999), 
and continuum, "hydrodynamic" versions of this model 
have been introduced (Toner and Tu, 1995, 1998; Topaz 
et al., 2006). Recently, Couzin et al. (2002, 2005) have 
introduced a model to provide insights into the mecha- 
nism of decision making in biological systems, which re- 
produces many important observations made in the field, 
and provides new insights into these phenomena. A re- 
view for various models can be found in Parrish et al. 
(2002) and Czirok and Vicsek (2001). 

The models of Couzin et al. (2005), and most other 
such models, incorporate various detailed mechanistic 
steps. These shed light on the role of leadership and 
imitation, and produce a number of surprising results, 
such as the influence that a few "informed" individuals 
can have on large collectives. What is needed now are 
efforts to simplify those models, and to show especially 
what properties of the microscopic simulators are essen- 
tial to explain that behavior. For some models, closure 
schemes are available (Flierl et al., 1999); but more gen- 
erally, though we may suspect that closures exist, we can- 
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not derive explicit expressions for them. In such circum- 
stances, we need methods such as those used in this pa- 
per; we perform the coarse-grained dynamical analysis by 
circumventing the derivation of governing equations, us- 
ing an equation-free computational approach (Thcodor- 
opoulos et al., 2000; Kevrekidis et al., 2003). A particular 
goal is to understand how much of the specific spatial de- 
tail is fundamental to the behavior. But turning to the 
Kuramoto-type approximation, where the interaction is 
assumed to be global, we deliberately ignore local effects. 
To the extent that the models fail to explain observed 
types of behavior, we will need to turn next to more de- 
tailed models. 

Most of previously proposed models concern popu- 
lations of homogeneous (or indistinguishable) individu- 
als. Furthermore, the dynamical analysis in the liter- 
ature is often limited to a small subset of the entire 
parameter space, and a systematic classification of pos- 
sible global dynamics remains elusive. In the current 
paper, we study the coarse-grained alignment dynamics 
of individual-based animal group models. The measure- 
ment of the mean angular deviation of fish schools (e.g. 
clupeids and scombroids; see Atz, 1953; Hunter 1966) 
showed that it varies continuously from no alignment 
to practically perfect alignment. We account for this 
continuous change by heterogeneity ("quenched noise"; 
characterized by parameters of random variables drawn 
from a prescribed distribution function) and the coupling 
strength. Our approach is flexible in that the heterogene- 
ity can be introduced in various places in the model, and 
the way we analyze different heterogeneity cases does not 
require any significant modification. 

The rest of the paper is organized as follows: Mod- 
els for homogeneous and heterogeneous animal groups 
are described in Sees. Ill Al and IIIBI and our approach, 
equation-free polynomial chaos, is explained in Sees. Ill CI 
and lll D1 Coarse-grained dynamical analysis and its com- 
parison with fine-scale dynamics, for a system of two in- 
formed individuals and a large number of heterogeneous 
uninformed individuals, are presented in Sec. IIIII The 
case of two groups of heterogeneous informed individu- 
als is presented in Sec. IIVI We conclude with a brief 
discussion in Sec. IVl 



II. MODELS AND METHODS 

A. A "minimal" model for identical individuals 

We briefly discuss a "minimal" model proposed by Na- 
bet et al. (2006), which we extend in our study. It con- 
cerns the alignment dynamics of a homogeneous popu- 
lation of indistinguishable N individuals with two sub- 
groups of informed individuals ( "leaders" ) with popula- 
tions iVi and A'2 respectively and -/V3 uninformed indi- 



viduals ("followers"), where iV = iVi + A^2 + A^a 
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Here ^pk characterizes the average direction of the indi- 
viduals in each of the two informed subgroups for A: = 1, 2 
and the average direction of the uninformed individuals 
for k — 3. Ofc is the corresponding informed, preferred 
direction {Qi can be set to zero without loss of general- 
ity) and K{> 0) is the coupling strength. This minimal 
model corresponds to the reduced system of the following 
system of A^ individuals (Nabet et al., 2006): 
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where the angle 9j characterizes the direction in which 
the jth individual is heading (we will refer to it as "orien- 
tation"). The average direction tpk is defined as the angle 
of the average of the phasors (when each individual's dy- 
namical state is considered as a phasor of unit radius and 
a phase angle) of the individuals in the fcth subgroup; pk 
is the magnitude of the average of the phasors. Formally, 
this is written as 
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The large population model in Eq. ^ has a separation of 
time scales. Individuals within each subgroup synchro- 
nize quickly, i.e., pk quickly converges to 1. The slow 
dynamics are described by the reduced system (Eq. ^) 
where the variables V'fc characterize the lumped behavior 
of each of the three subgroups. 

It is assumed that the alignment (orientational) dy- 
namics are independent of the translational counterpart 
(Sepulchre et al., 2005); hence, the dynamical state of 
an individual can be characterized by its orientation. 
The functional form for mutual interaction is borrowed 
from the well-known Kuramoto model (Kuramoto, 1984), 
a prototypical model for coupled nonlinear oscillators. 
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This simplified global interaction model is consistent with 
an observation that the strongest correlations are ob- 
served between the (speed and) direction of the individ- 
ual and the average (speed and) direction of the entire 
school (Partridge, 1982): In the mean-field form of the 
Kuramoto model, the interaction term can be rewritten 



1 ^ 
N ^ 



smit/; 



,) = rsin(V' - 6'j), 



(3) 
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In this alternate expression, the dependence on the dif- 
ference between the individual direction and the average 
direction stands out explicitly. In the absence of coupling 
[K = 0), each leader eventually heads for its preferred 
direction. Nontrivial dynamical behavior for the minimal 
model (Eq. ^) are studied in Nabet et al. (2006); bi- 
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furcations are analyzed for the global phase space in the 
case iVi = N2 and N3 = 0. 



B. Extension to heterogeneous populations 

The aforementioned models concern populations of 
homogeneous subgroups, where the individuals in each 
subgroup quickly synchronize, nearly perfectly (p^ ~ 1), 
during the initial transients (Nabet et al., 2006). In 
the more general case, the mean angular deviation of 
fish schools is finite (Atz, 1953; Hunter, 1966), which is 
not captured in this "minimal" model. We extend the 
model to account for the distribution of directions within 
schools, assuming it arises from the heterogeneity among 
the group members. We introduce the heterogeneity in 
the following two ways: 

(I) Two leaders and many heterogeneous followers — 
First we consider the cases when the population consist of 
two leaders (which possibly represent lumped behavior of 
groups of homogeneous leaders) and N 1) followers: 
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for l<i< N, 



where the heterogeneity is accounted for through the 
tendency to deviate from the average direction, char- 
acterized by oji, an i.i.d. random variable drawn from 
a prescribed distribution function g{uj) (of standard 
deviation cr„ with mean value zero). For notational con- 
venience, we drop a subscript of a variable to represent 
a random variable of a proper length (c/. oji and lo). 
As Qi can be set to zero without loss of generality, 82 
and K are control parameters. In the current study, we 
consider g{uj) to be Gaussian, but our analysis is not 



limited to this particular choice. 

(II) Two groups of heterogeneous leaders — Secondly, 
we consider two groups of heterogeneous leaders without 
any followers, focusing only on the dynamics among lead- 
ers. The heterogeneity is accounted for by introducing 
randomness in the angles preferred by the leaders. The 
orientations of the leaders in each group are denoted by 
Xi's and (j)iS, (of sizes A^i and A^2) respectively: 
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for 1 < i < A^i, 
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FIG. 1: Direct integration of a system of two leaders (open 
circles; dashed lines indicate preferred angles) and 300 follow- 
ers (dots), initialized from uniformly distributed orientations 
with randomly assigned heterogeneity variable (i.e., no initial 
correlations between 6 and ui), is shown for an initial transient 
[(a) to (c)], and for much longer time scales [(d) to (f)]. Insets 
illustrate time evolution of the followers' orientations on the 
9 — u> plane, where strong correlations develop during a short 
time t ~ 10. After that, the leaders and followers, the latter 
effectively as a "unit" , slowly drift to the stable steady state. 
It takes of the order of t ~ 10^ for the system to asymptoti- 
cally converges to this final state. {K — 1.0; Q2 = 7r/4). 



where the preferred angles Xi and <I>i are randomly drawn 
from prescribed distributions gi{X) and (72 (^') (i-e., i.i.d. 
random variables of standard deviations ax and cr$), re- 
spectively. We set < A:" 0, and will vary K {> 0) and 
< $ > (G [OjTt]) as control parameters (and investigate 
some cases of different values of ax and a<s> in Sec. lIVB]) . 

C. Wiener polynomial chaos 

The Kuramoto model, a paradigm for all-to-all, phase- 
coupled oscillator models, has been extensively studied 
and used to shed light on many synchronization phenom- 
ena (Kuramoto, 1984; Acebron et al., 2005, and refer- 
ences therein). This model has the property that, in the 
full synchronization regime (of large enough K values), 
phase angles become quickly correlated with (or "sorted" 
according to) the natural frequencies during the initial 
short transients (Moon et al., 2006). Similar correlations 
(between the angles and heterogeneity random variables) 
are expected to arise in the current model (which is in- 
deed the case, as will be shown later in Fig. [T|) , since the 
coupling term is qualitatively similar. As in Moon et al. 
(2006), we choose expansion coefficients in Wiener poly- 
nomial chaos as coarse-grained "observables" , to explore 
low-dimensional, coarse-grained dynamics. 

Wiener (-Hermite) polynomial chaos was introduced by 
Wiener (1938), who represented a random process in 
terms of functional expansions of Wiener process (his- 
torically, this method has been termed as polynomial 



"chaos", because of its initial usage on homogeneous 
chaos, such as turbulence and Brownian motion, rather 
than the nature of the method). Ghanem and Spanos 
(1991) later extended this idea to treat random processes 
as functional expansions of random variables, or elements 
in the Hilbert space of random functions, in which a spec- 
tral representation in terms of polynomial chaos is iden- 
tified. The projections (or coefficients) on the polyno- 
mial base then can be determined through a Galerkin 
approach. This method was subsequently applied in 
uncertainty quantification of various problems (e.g., see 
Ghanem, 1999; Jardak et al., 2002), and has been ex- 
tended to general situations using the Askey scheme (Xiu 
et al., 2002; now known as generalized polynomial chaos). 

In this method, dependent random variables {9 of the 
followers for the case (I), and x a-nd (p for the case (II)) 
are expanded in polynomials of independent random 
variables (w, or X and $) using appropriately chosen 
basis functions. Details for the two cases are as follows: 

(I) Two leaders and many heterogeneous followers — 
For convenience, we introduce the unit Gaussian random 
variable ^ = Loja^. Using this newly defined variable, 
we expand 9(uj,t) (i.e. 9{S,,t)) in Hermite polynomials 
of C [HoiO - 1, H,iO = ^, H^iO =e-h HsiO = 

p 

0(^,t) = ^a„(i)i7„(C), (6) 

n=0 

where p is the highest order retained in the truncated se- 
ries, Hn is the nth Hermite polynomial, and the a^s are 
the expansion coefficients which will be referred to simply 
as "chaos coefficients" in this paper. Wiener polynomial 
chaos, utilizing Hermite polynomials as basis functions, 
is the appropriate choice for Gaussian random variables 
that we consider in the present study. The probability 
density function of the Gaussian random variables ap- 
pears as the weighting function of Hermite polynomials, 
and the Hermite polynomial expansion is suggested to 
converge exponentially for Gaussian processes (Lucor et 
al., 2001). For other random variables, use of different 
basis functions (for instance, Legendre polynomials for 
uniform random variables) has been suggested for fast 
convergence, which is the basis of the development of the 
generalized polynomial chaos (Xiu et al., 2002). 

We choose the first few nonvanishing chaos coefficients 
a„'s, as well as the orientations of the leaders (■0i and 
-02), to be the coarse-grained "observables". Due to 
symmetry, all the even order a^'s vanish, except for the 
zeroth order ao that corresponds to the average direction 
of the followers. Geometrically, ai and 013 respectively 
represent a measure for the linear order spread of the 
angles (the "slope" between 9 and oj) and the cubic 
order measure. In the continuum limit [N — > 00), the 
chaos coefficients can be exactly determined using the 
orthogonality relations for Hermite polynomials. How- 
ever, in the finite cases of single realization we consider. 
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N ~ O(IO^), those relations hold only approximately, 
and the coefficients are evaluated using least squares 
fitting, following Moon et al. (2006). 

(II) Two groups of heterogeneous leaders — In the sec- 
ond case, we expand x in terms of X and $, re- 
spectively: 

p 

X = ^a„i7„(C), 

n=0 
P 

^ = Y^f^nHniv), (7) 

n=0 

where the chaos coefficients a and /3 are the coarse "ob- 
servables" of our choice, iJ„'s are Hermite polynomials 
(for Gaussian gi and 52), and C = X/a^ and r/ = ^/cr^ 
are unit Gaussian random variables. 



D. "Equation- free" computational approach 

A prerequisite to coarse-grained dynamical analysis 
(which is the main goal of the current study) is, in a tradi- 
tional sense, an explicit derivation of coarse-grained gov- 
erning equations. In principle, such equations for chaos 
coefficients, in the continuum limit (N —> 00), might be 
obtained through a stochastic Galerkin method (Ghanem 
and Spanos, 1991). 

In the present study, we do not even attempt to derive 
such equations. We circumvent their derivation by us- 
ing an equation-free multiscale computational approach 
(Theodoropoulos et al., 2000; Kevrekidis et al., 2003, 
2004). This approach enables us to explore the coarse- 
grained dynamics without the assumption of the contin- 
uum limit. The premise of this approach is that coarse- 
grained governing equations conceptually exist, but are 
not explicitly available in closed form. This approach is 
based on the recognition that short bursts of appropri- 
ately initialized microscopic (fine-scale) simulations dur- 
ing a time horizon AT and the projection of the re- 
sults onto coarse-grained variables, say x, result in time- 
steppers (mappings) for those variables $at (which is 
effectively the same as the discretization of unavailable 
equations): 

Xn+l = ^AT{Xn)- (8) 

One then processes the results of the short simulations to 
estimate various coarse-grained quantities (such as time 
derivatives, action of Jacobians, residuals) to perform rel- 
evant coarse-grained level numerical computations, as if 
those quantities were obtained from coarse-grained gov- 
erning equations. For instance, one can integrate unavail- 
able governing equations in time (coarse projective inte- 
gration; see below), or compute the steady states of the 
above coarse time-stepper, by utilizing fixed point algo- 
rithms (such as Newton- Raphson or Newton-GMRES). 



Equation-free computations consist of the following 
steps: 

1. Identify coarse-grained variables ("coarse observ- 
ables") that sufficiently describe both the micro- 
and macroscopic dynamics; in our study, they are 
a„'s (and /3n's). For convenience, we denote the 
microscopic (macroscopic) descriptions by 9 (a). 

2. Choose an appropriate lifting operator /j.^, which 
maps a to one (or more) consistent description(s) 
9 (for the purposes of variance reduction and 
ensemble-averaging). In the current study, this can 
be achieved by using the relations in Eqs. ([6]) and 
([7|: once random variables are drawn, these rela- 
tions are used to obtain corresponding 9. 

3. Starting from lifted initial condition(s) 9{to) = 
/iL(a(io)), run the microscopic simulator to obtain 
9{to + AT) at a later time (AT > 0). 

4. Use an appropriate restriction operator Ada (least 
squares fitting, in the current study) which maps 
the microscopic state(s) to the macroscopic descrip- 
tion a{to + AT) = MR{9{to + AT)), which effec- 
tively results in time series of coarse observables, 
or their coarse time-stepper $at; «(^o + AT) = 
*AT(a(io))- 

5. Apply desired numerical techniques using the 
coarse-grained variables obtained from the step 4. 
and repeat some of the above steps as needed. 

An extensive discussion can be found in Kevrekidis et al. 
(2003, 2004). 



III. RESULTS FOR CASE I 

Direct integration of the "fine-scale" model of Eq. Q 
in the strong coupling regime {K — 1.0, a^^ = 0.1), 
started from randomly assigned orientations and the het- 
erogeneity variable (the latter is a Gaussian random vari- 
able) , illustrates that a strong correlation between 9 and 
u) develops during a short, initial transient time; the ori- 
entations of the followers quickly become a monotoni- 
cally increasing function of their heterogeneity variable 
(Fig. [T]), after which they slowly drift as a "unit" until 
they settle down in the final steady state. During the lat- 
ter slow drift, the system can be described as two leaders 
and a single "clump" of followers, whose coarse-grained 
states can be successfully described by a small number of 
chaos coefficients. A similar time scale separation exists 
in the model of homogeneous populations. In this case, 
followers quickly collapse asymptotically to the same di- 
rection (Nabet et al., 2006). 




FIG. 2: (Color online) Accelerated computation of sta- 
ble steady states via coarse projective integration using five 
coarse-grained variables, shown here for two different time 
scales {K — 1.0; O2 ~ 7r/4). Initially all the values are as- 
signed to be 0. Both ai and 0:3 reach their steady state values 
relatively quickly (see (a)), while the others are slowly vary- 
ing (see (b); they are still varying nt t = 500). Dots represent 
the time intervals during which short direct integration is per- 
formed (and restricted), in the course of the projective inte- 
gration using forward Euler method. Solid lines represent the 
trajectories of direct full integration during the entire time. 
Higher efficiency can be achieved by optimally choosing the 
time horizon for the direct integration, the projection stepsize, 
and projection method. 



A. Computations of steady states 

We begin by accelerating the approach to a stable 
steady state using an equation-free algorithm, the coarse 
projective integration method (Gear and Kevrekidis, 
2003). In contrast to a conventional, direct integra- 
tion of the full fine-scale model during the entire time 
(until sufficient convergence to stable, final states), this 
method exploits smoothness in the coarse variables (esti- 
mated through a direct integration during a short time) , 
in order to extrapolate and take a large projective time- 
step (compared to the original integration time-step size). 
This saves computational effort. The procedure consists 
of (z) lifting (appropriate initialization of the fine-scale 
simulator, an integrator of Eq. consistent with pre- 
scribed coarse-grained values), (ii) direct integration of 
the microscopic simulator during a relatively short time 
interval (but long enough to accurately estimate local 
coarse-grained time derivatives), (Hi) restriction (of fine- 
scale description onto coarse-grained variables), and (iv) 
taking a projective step (using a traditional numerical in- 



FIG. 3: (Left panel) A bifurcation diagram observed on one 
leader V2, computed using AUT02QQQ {K = 1.0). Solid 
(dashed) lines represent stable (unstable) branches. There 
exist a few other unstable branches that are not shown here. 
At some critical value of O2, an unstable state in the up- 
per branch undergoes a forward pitchfork bifurcation; two 
unstable states coincide. The lower branch of "trivial" solu- 
tions does not exhibit any bifurcation. (Right panels) Snap- 
shots of two (symmetric) stable states in the bistable regime 
(G2 = 2.0), marked by dots in the left panel. 



tegration scheme such as forward Euler). The computa- 
tional payoff of this method depends on the ratio between 
a short direct integration time interval, the projective 
time-step size, and the computational effort required for 
lifting/restriction procedures (see e.g., Rico-Martinez et 
al., 2004). More importantly, successful computation of 
steady states through this method naturally attests to 
the validity of the chosen coarse-grained observables in 
describing both fine-scale and coarse-grained states. 

Projective integration using five coarse-grained vari- 
ables {ipi,ip2, and the first three non-vanishing a„'s; ao, 
ai, and a^) follows virtually the same trajectories of the 
full, direct integration (Fig. ^ , even if to is newly drawn 
at each lifting; the agreement is even better if the same 
uj were used (hence the dynamics are fully determinis- 
tic). Both lifting (simply using Eq. ^) and restriction 
(least squares fitting) operations require minimal compu- 
tational efforts. Therefore, the computational efficiency 
in the present case is nearly exclusively determined by 
the projective step size, which is a factor of about four in 
Fig. [51 with a more sophisticated projection algorithm, a 
higher efficiency can be obtained. We see that both ai 
and Q!3 reach their steady state values quickly {t ~ 5), 
showing that the correlation between 9 and uj are fully de- 
veloped by then. However, the other chaos coefficient aQ 
(representing the average direction) slowly drifts towards 
the steady state, and so do ipi and ip2 (note that it is still 
varying at t = 500); the computation of an asymptotic, 
steady state requires a very long time integration. 

Direct integration (including projective integration) 
cannot compute unstable steady states and are inappro- 
priate for stability computations and parametric bifur- 
cation studies. Both stable and unstable steady state 
values can be systematically (and much more efficiently 
than the projective integrations) computed by applying 
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TABLE L A coarse steady state computation ai K — 1.0 and 02 — it/4: for A'' — 300, using the Newton-GMRES method. 
Values at each iteration have been averaged over an ensemble of 100 reahzations. The last column shows relative residuals. 
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coarse-grained fixed point algorithms to the steady state 
condition of Eq. i.e., x — ^at{x) = 0, in much 

lower dimension than that of individual level. We use 
the coarse Nevirton-GMRES (Kelley, 1995), a matrix- free, 
method to compute coarse-grained fixed points. We ob- 
serve that the algorithm accurately converges within a 
few steps (Tab. |l]); the converged values are accurately 
consistent with the restricted values of the fixed point so- 
lution of the detailed (i.e., (A^-f 2)-dimensional) problem, 
within prescribed convergence tolerance. By combining a 
coarse fixed point algorithm with pseudo-arclcngth con- 
tinuation (Keller, 1987), we numerically compute coarse- 
grained bifurcation diagrams in the following sections. 
The computational efficiency of the coarse fixed point al- 
gorithm varies with the choice of the initial guess for the 
iteration. With a totally uneducated guess, it could take 
even longer than the direct integration (note that the lat- 
ter never computes the exact solutions); however, during 
the continuation computation shown below, a good ini- 
tial guess is always available from the previous parameter 
value (s), and even several orders of magnitude of compu- 
tational efficiency can be achieved. 



B. Types of fine-scale dynamical behavior 

We first analyze the detailed ( A^ -I- 2)-dimensional fine- 
scale model in the full synchronization regime, in or- 
der to obtain insights on fine-scale dynamics to be com- 
pared with our coarse-grained analysis below. We use 
AUTO2000 (Doedel et al., 2000) to compute the fine- 
scale bifurcation diagrams as functions of O2 at a fixed 
value of K; only projections for one leader {^2) are shown 
in Fig. [3] and for one follower in Fig.|4](a). All the other 
followers exhibit essentially the same dynamical behav- 
ior as the one shown here (except for some quantitative 
differences) . 

The interaction between the individuals causes the 
steady state directions of the leaders to deviate from 
the preferred angles and 62, respectively. Such devia- 
tion can occur in two directions, either toward the region 
bounded by [0, 62] (an "obvious" steady state where fol- 
lowers are directed in between the directions of the lead- 
ers; see Fig. [3] (b)) or the other way around (e.g., Fig. [3] 



(a) 




FIG. 4: (a) A bifurcation diagram observed on one (arbitrar- 
ily chosen) follower, as a function of 02 {K = 1.0), computed 
using AUTO2000 (the same case as in Fig. [3|. Superscript 
'f has been added to emphasize that this is the orientation 
of a follower. A few other existing unstable branches are 
not included here. The upper branch undergoes a pitchfork 
bifurcation and becomes stable, (b) A coarse bifurcation di- 
agram observed on qq (average direction), obtained by the 
coarse Newton-GMRES method with pseudo arc-length con- 
tinuation. Only a blowup around the bifurcation point is 
shown. Coarse-grained dynamics exhibit the same structure 
as in the fine-scale level. Filled (open) circles represent stable 
(unstable) steady states. 



(a)). The analysis shows that for small 82 values only 
the former state is stable, while for large values, both of 
these states become stable. The branches for "obvious" 
stable steady states, which correspond to lower straight 
solid Hnes, exhibit no bifurcations (see Figs.[3]and|4](a)). 
On the other branches, forward pitchfork bifurcations 
at some critical value of O2 give birth to another sta- 
ble branch (a state on this stable branch is shown in 
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Fig. [3] (a)), as well as two unstable asymmetric solution 
branches, hence the population becomes bistable. The 
critical value of Q2 for the onset of the bistability de- 
pends on K (precisely speaking, K/aui); the critical value 
is 82 0.45 (2.2) at K ^ 1.0 (0.5). As K decreases fur- 
ther, the critical value monotonically increases until fully 
synchronized steady states lose stability at some critical 
value of K. 



C. Coarse-grained dynamics 

We now compute coarse-grained steady state solutions. 
A coarse-grained bifurcation diagram for ao (represent- 
ing the average direction of the followers) is compared 
with the corresponding diagram observed for one fol- 
lower, in Figs, m (b) and (a); (b) is a blowup of the 
region around the bifurcation. Both sides of the bifur- 
cation point can be described by the same set of coarse- 
grained observables, which clearly summarize group level 
dynamical behavior of the followers before and after the 
bifurcation. 

As K decreases in the Kuramoto model, oscillators get 
desynchronized (Kuramoto, 1984), starting with the os- 
cillator with the maximum value of \LOi\ (the "extreme" 
oscillator) through a saddle-node (actually a "sniper") bi- 
furcation on a limit cycle (Moon et al., 2006). We expect 
the same type of bifurcation to occur in this model. How- 
ever, when we try to compute the coarse-grained steady 
states as functions of K using the previously mentioned 
five coarse variables (via coarse Newton-GMRES method 
and pseudo arc-length continuation, neither a bifurcation 
nor an unstable branch is appropriately identified. The 
computation, initialized at large K steady states, accu- 
rately follows stable branches down to some critical value 
of K (where the transition occurs) , and then fails to con- 
verge. Our coarse-grained observables are not sufficient 
to describe the states on the "other side" of the bifurca- 
tion point, as we will explain below. 

A fine-scale bifurcation diagram (computed using 
AUTO2000) obtained by starting from a stable steady 
state on the lower branch in Fig. [3] is shown in Fig. [5] 
(a). Here the diagrams for two leaders and only a few 
followers, including the extreme one, are shown. We find 
that both stable and unstable branches for each angle 
nearly coincide for all the individuals (see inset of Fig. [5] 
(a)), except for the extreme one. As the difference be- 
tween stable and unstable branches (at the same value 
of K) is appreciable only when observed on this extreme 
oscillator, a smooth mapping between 6 and uj does not 
prevail for unstable states, and the previously used chaos 
coefficients are not appropriate any more. 

Taking these observations into account, it is easy to 
remedy the situation as follows: The fact that stable and 
unstable branches nearly coincide, discounting the ex- 
treme follower, suggests that all the individuals except 
for the extreme follower can be again described by the 
same set of chaos coefficients. Thus we treat the orienta- 
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FIG. 5: (Color online) (a) A bifurcation diagram observed on 
a few followers, including the one with the maximum of \uji\ 
(the "extreme" follower), as a function of K (Q2 = 7i"/4), com- 
puted using AUTO2000. A critical value where the extreme 
individual gets desynchronized corresponds to a saddle-node 
bifurcation point on a limit cycle (a "sniper" bifurcation). Ex- 
cept for the extreme follower, stable and unstable branches 
nearly coincide (see the inset), (b) In order to capture the 
fine-scale bifurcation, the angle of the extreme follower has 
to be discounted from the chaos expansion and considered as 
an extra coarse-grained variable (see text). We distinguish 
these chaos coefficients (from the ones used so far) by adding 
a prime. It was computed via the coarse Newton-GMRES 
method with continuation. 



tion of the extreme one separately (introducing it as an 
additional coarse-grained variable) , and discount it from 
the polynomial chaos expansion. (From the fact that 
the extreme follower gets desynchronized at the transi- 
tion, one can also intuitively see that followers have to 
be considered as a combination of a clump of synchro- 
nized "bulk" and a separate, extreme one.) We compute 
the solutions with continuation, using this new set of six 
coarse variables, which captures the bifurcation and ap- 
propriately describes the unstable steady states (Fig. [5] 
(b)); we have analyzed exactly the same realization used 
in Fig. [5] (a) for direct comparison. When bifurcation 
diagrams are computed for ensembles of many realiza- 
tions, an uncertainty will arise in the exact quantification 
of the bifurcation point, due to the fluctuation of finite- 
dimensional random variables among realizations, while 
the results are qualitatively the same as those of a single 
realization (Xiu et al., 2005). 

The coarse bifurcation results shown in Fig. [5] (b) il- 
lustrate that the steady state directions of the leaders 
and the average direction of followers (ctg, discounting 
the extreme one; a prime is added to distinguish it from 
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FIG. 6: Lines: A bifurcation diagram observed on tpi in 
the minimal model (the first two ODEs in Eq. ([ij with 
A'^i = A^2 = 1, = 0), for varying O2 at a fixed value 
oi K = 2.4, obtained by AUTO2000. For large enough pre- 
ferred angles (02/7r >~ 0.7), the system becomes bistable 
through a forward pitchfork bifurcation. Circles: A coarse 
bifurcation diagram near the pitchfork bifurcation, observed 
on Qfo, as a function of < $ > /tt, computed by the coarse 
Newton-GMRES method. Filled (open) circles represent sta- 
ble (unstable) steady states. 



the previously used notation) are virtually the same for 
a range of K. Only higher order chaos coefScients (only 
a'l is shown in Fig. [5]) appreciably vary as a function of 
K, which means that individuals spread more widely as 
K decreases, until the extreme one eventually starts to 
oscillate freely, while the average steady state direction 
remains the same. 



IV. RESULTS FOR CASE II 
A. Dynamics of statistically similar groups 

Here we explore both the fine-scale and coarse-grained 
dynamics of a model for two groups of heterogeneous lead- 
ers (with no followers) shown in Eq. ([5]), and compare the 
results of the two different scales. One notable difference 
from the Kuramoto model is that "oscillators" in Eq. ^ 
do not have finite random variables (natural frequencies), 
hence there is no onset of the synchronization that occurs 
at a finite value of K (or, they can be alternatively seen 
as Kuramoto- like oscillators of zero natural frequencies, 
which result in the onset at K — 0, hence they get syn- 
chronized for all K values). The analysis of the minimal 
model (the first two of Eq. ^ with Ni ^ N2, N3 ^ 0) 
reveals that for large enough Q2 (>~ 7''/2) the system 
exhibits bistability for a certain range of K (Nabet et 
al., 2006), as in the previous case in Sec. IIIII Here we 
will vary < $ > as the main parameter for two differ- 
ent values of K. For large coupling strengths {K > 2.0), 
the bistability in the minimal model appears through a 
forward pitchfork bifurcation, when &2 is varied as a pa- 
rameter (Fig. [5]) . This minimal model can be seen as a 
special case of the current model, where both X and <i> 
are assumed to be delta functions and each group consists 
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FIG. 7: Coarse-grained bifurcation diagrams observed on the 
first two chaos coefficients: (a) ao; the average direction of the 
first group of leaders, and (b) ai ; the "slope" between x BJad 
X, as functions of < $ >. These are blowups of the region 
around the forward pitchfork bifurcation point in Fig. |6] 



of identical individuals. 

We begin by asking whether our model for hetero- 
geneous groups exhibits similar types of dynamical be- 
havior. One can also do accelerated computations of 
steady states using the coarse projective integration, but 
here we skip such computations and present only the 
coarse bifurcation analysis results. Coarse bifurcation 
diagrams obtained through the coarse Newton-GMRES 
method (Kelley, 1995) and pseudo arc-length continua- 
tion (Keller, 1987) (for Gaussian distributions of X and 
$; CTA- = (T$ = 0.1, Ni = N2 = 100) show that the het- 
erogeneous groups indeed exhibit the same qualitative 
type of coarse dynamical behavior around the pitchfork 
bifurcation point (Fig. [7]). As we consider symmetric 
unimodal distribution functions, all the even order chaos 
coefficients (except for ao and (3q) virtually vanish. The 
diagram for ao of the first group < x > (average di- 
rection) exhibits reasonably good quantitative agreement 
with the corresponding diagram for the minimal model, 
within fluctuations of finite-size random variables, shown 
in Figs. [5] and [5] It is interesting to note that at the criti- 
cal point, all the followers are headed for the same direc- 
tion (ai — 0, which corresponds to the "slope" between 
X and x); see Fig. [7](b). 

The Hermite polynomial expansion converges so 
quickly that the expansions can be accurate even when 
truncated at the third order. Due to the reflection sym- 
metry (about < $ >/2), /3 coefficients have similar struc- 
tures as the a ones, after proper reflection and transla- 
tion. Only results on a are presented here. As the cou- 
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FIG. 8: Lines: A bifurcation diagram for the minimal model 
of two leaders, for varying O2 at a fixed value of if = 1.8, 
obtained by AUTO2000. For large enough preferred angles 
(02/7r >~ 0.75), the system becomes bistable, but the nature 
of the bifurcation is different from that of higher K value 
cases (a saddle-node vs. a pitchfork bifurcation; see Fig. [S]). 
Circles: A coarse bifurcation diagram observed on the average 
direction of the first group of leaders {010) around the saddle- 
node bifurcation, as a function of < $ > /tt, computed via 
the coarse Newton-GMRES method with continuation. Filled 
(open) circles represent stable (unstable) steady states. 



pling strength decreases across K — 2.0, the nature of the 
bifurcation changes (from a pitchfork) to a saddle-node 
bifurcation (Fig. ^ at K — 2.0, which also occurs in the 
model for homogeneous populations; the nature of the 
transition between these different bifurcations, a higher 
codimcnsion bifurcation, has been discussed in Nabet et 
al. (2006). 
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FIG. 9: (Color online) Coarse-grained bifurcation diagrams 
near a turning point in Fig. [S] for X distributions of three 
different widths (the standard deviation ax = 0.1 for circles; 
0.2 for squares; 0.3 for triangles), obtained via the coarse 
Newton-GMRES method and continuation. The standard 
deviation for the second group, a-s>, is kept the same at 0.1 
{K = 1.8, iVi = iV2 = 100). Filled and open symbols rep- 
resent stable and unstable states, respectively, (a) The first 
chaos coefficients ao (average direction of the first group) are 
nearly the same for the three cases. The difference between 
the cases becomes apparent in higher order coefficients that 
reflect the degree of spreading; see Qi in (b). 



B. Statistically different groups 

So far we have considered statistically similar groups, 
namely iVi — N2 and ax — 0"$; they differed only by 
average preferred directions. It is natural to ask how 
the dynamics change as the parameters concerned with 
the distributions (for the preferred directions) are var- 
ied. It is readily expected that the essential dynamics 
of two different-size groups can be reflected in the mini- 
mal model using two different coupling strengths, which 
is considered in Nabet et al. (2006). Here we consider 
only the cases with varying width of the distributions 
(cr$ ^ (^x), which has no analog in the minimal model. 

Coarse bifurcation diagrams for three different Gaus- 
sian distributions for X (ax is varied while cr$ is kept at 
0.1; see Fig. ^ show that the average directions (ao's) 
hardly vary with the width of the distributions; the pri- 
mary parameter that affects on the average direction is 
the group size. For the distributions of different widths, 
the fixed point computation with continuation fails to 
converge at different values of ao's; points marked by ar- 
rows m Fig. M are the last points the Newton-GMRES 
computations converged in each case, when approached 
from the stable branches. Such a failure of conver- 
gence can be expected, because the steady states on this 



unstable branch overlap with another nearby unstable 
branch (which is not shown in this figure, but was shown 
in Fig. [5]); characterizing the distribution with a few 
Wiener chaos coefficients does not provide an accurate 
description any more. The differences between the three 
cases (of different distribution widths) manifest them- 
selves clearly in higher order chaos coefficients. While 
the average behavior remains nearly the same (Fig. [S] 
(a)), individuals in the group spread more widely (as re- 
flected in ai and higher order coefficients; Fig. [21(b)), as 
the width of the random variable (distribution) increases. 



V. CONCLUSIONS 

We have demonstrated a computational venue (an 
equation-free polynomial chaos approach) to study 
coarse-grained dynamics of individual-based models ac- 
counting for the heterogeneity among the individuals 
in animal group alignment models. We considered fi- 
nite populations of (I) two "leaders" (which have di- 
rect knowledge on preferred directions) and N{'^ 1) 
uninformed, heterogeneous "followers" , and (II) two 
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groups of heterogeneous "leaders". We explored the 
coarse-grained, group level (low-dimensional) dynamics 
using the polynomial chaos expansion coefhcients as 
coarse-grained observables; these observables account for 
rapidly developing correlations between random vari- 
ables, and sufficiently specify both fine-scale and coarse- 
grained (group-level) dynamical states. 

All the analysis in our study was done expressively 
avoiding the derivation of coarse-grained governing equa- 
tions, following a nonintrusive, equation- free computa- 
tional approach wrapped around the direct system simu- 
lator. It should be noted that we have not assumed that 
N is infinitely large (so-called the "continuum limit"). 
Our approach can be used for systems of any finite, large 
number of populations, and it can be equally applied 
to various types of random variables (following general- 
ized polynomial chaos) and/or various heterogeneity. We 
compared our results with those of minimal models that 
do not account for heterogeneity among the individuals. 
They show good agreement in the lowest order (i.e., av- 
erage directions) , which clearly highlights the correspon- 
dence between the individual- and group-level dynamics 
(Figs. [H] and [5]). Indeed this implies that the results in 
Nabet et al. (2006), where no heterogeneity is explic- 
itly accounted for, are more robust than demonstrated 
in that paper alone. 

In order to analyze different coarse-grained bifurca- 
tions, it became necessary to use different sets of coarse- 



grained variables, even if the model is the same in the 
fine-scale level (Fig. [5]). This clearly shows that an ap- 
propriate choice of coarse-grained observables (in terms 
of which one can obtain useful closures) is an essential 
step; different coarse-grained observables are required, 
as the same fine-scale model closes differently. 

In the present study, we assumed that the orientational 
dynamics can be separated from their translational coun- 
terpart, and considered the simplest nontrivial cases of 
all-to-all ( "all- visible" ) , sinusoidal coupling. Our future 
work will involve the incorporation of the translational 
dynamics and more complicated coupling/network topol- 
ogy, including heterogeneous couplings. Our work pre- 
sented here is the first step of our effort toward the de- 
velopment of more detailed (and biologically more plau- 
sible) models and their coarse-graining. 
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